Análisis geospacial con Python

Introducción

Paquetes de análisis espacial:

  • Geopandas: paquete para datos vectoriales (equivalente a sf)

  • Xarray: paquete para datos multidimensionales (equivalente a stars?)

  • Rasterio: interfaz a GDAL para datos raster (equivalente a terra, pero de bajo nivel)

  • Rioxarray: puente entre Xarray y Rasterio.

Paquetes de visualización interactiva:

  • Leafmap: interfaz sencilla y unificada para mapas web de múltiples backends (folium, ipyleaflet, maplibre).

  • MapLibre: mapas 2D y 3D avanzados.

  • HyperCoast: visualización 3D de datos hiperespectrales con pocas líneas de código.

Paquetes de análisis especializado:

  • WhiteboxTools: herramientas de geoprocesamiento de análisis hidrológico, análisis del terreno, y análisis de datos LiDAR.

  • Geemap: puente a Google Earth Engine.

  • DuckDB: análisis de datos espaciales en base de datos de alto rendimiento.

  • Apache Sedona: computación distribuida para procesado de datos geoespaciales.

Desarrollo de aplicaciones:

  • Voila y Solara: creación de aplicaciones con Python y Jupyter.

  • Shiny: creación de aplicaciones y dashboards.

Vamos a comprobar que los paquetes se cargan correctamente y a crear un primer mapa:

## cargar paquetes
import geopandas as gpd
import rasterio
import xarray as xr
import rioxarray
import leafmap
import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
## crear mapa
map = leafmap.Map(center = [40, -100], zoom = 5, height = '500px')

## añadir mapas base
map.add_basemap('OpenTopoMap')
map.add_basemap('USGS.Imagery')

## mostrar mapa
map

Análisis vectorial con GeoPandas

Vamos aprender:

  • GeoDataFrame y GeoSeries

  • Leer/Exportar datos vectoriales

  • Operaciones comunes

  • Visualizar datos

  • Trabajar con CRS

Geopandas

Una diferencia con sf, es que un GeoDataFrame puede tener varias geometrías, pero solamente una está activa en cada momento. La geometría se accede a través del atributo .geometry. Vamos a crear un GeoDataFrame:

## datos
data = {
    'City': ['Tokyo', 'New York', 'London', 'Paris'],
    "Latitude": [35.6895, 40.7128, 51.5074, 48.8566],
    "Longitude": [139.6917, -74.0060, -0.1278, 2.3522]
}

## convertir a DataFrame
df = pd.DataFrame(data)

## convertir a GeoDataFrame
gdf = gpd.GeoDataFrame(df, geometry = gpd.points_from_xy(df.Longitude, df.Latitude))

## imprimir
print(gdf)
       City  Latitude  Longitude                  geometry
0     Tokyo   35.6895   139.6917  POINT (139.6917 35.6895)
1  New York   40.7128   -74.0060   POINT (-74.006 40.7128)
2    London   51.5074    -0.1278   POINT (-0.1278 51.5074)
3     Paris   48.8566     2.3522    POINT (2.3522 48.8566)

Leer datos:

## url a geojson
url = "https://github.com/opengeos/datasets/releases/download/vector/nybb.geojson"

## leer
data_geojson = gpd.read_file(url)

Exportar datos:

gdf.to_file('00_data\\test-data.geojson')
gdf.to_file('00_data\\test-data.gpkg', layer = "layer1")

gpd.list_layers('00_data\\test-data.gpkg')
name geometry_type
0 layer1 Point

Sistemas de Referencia de Coordenadas

Comprobar CRS:

print(data_geojson.crs)
EPSG:2263

Reproyectar o transformar CRS:

gdf_4326 = data_geojson.to_crs('EPSG:4326')

Medidas

Para calular medidas como áreas y distancias, necesitamos un CRS que esté en unidades como metros (p. ej. EPSG:3857).

## reproyectar
gdf_3857 = gdf_4326.to_crs('EPSG:3857')

## crear índice
gdf_3857.set_index('BoroName', inplace = True)

Calcular área:

## área
gdf_3857["area"] = gdf_3857.area

## perímetro en km
gdf_3857["perim"] = gdf_3857.length / 1000

Extraer características geométricas

## centroide
gdf_3857['centroid'] = gdf_3857.centroid

## bordes
gdf_3857['boundary'] = gdf_3857.boundary

Distancias: vamos a ver la distancia desde el centroide de Manhattan al resto de centroides

## centroide de referencia
manhattan_centroid = gdf_3857.loc['Manhattan', 'centroid']

## calcular distancia
gdf_3857['dist_to_manhattan'] = gdf_3857['centroid'].distance(manhattan_centroid)

Crear buffers:

## crear buffer (metros)
gdf_3857['buffer'] = gdf_3857.buffer(3000)

## visualizar
fig, ax = plt.subplots(figsize=(10, 6))

gdf_3857["buffer"].plot(
    ax        = ax,
    alpha     = 0.3,
    color     = "orange",
    edgecolor = "red",
    linewidth = 1,
    label     = "3km Buffer Zone",
)

gdf_3857.plot(
    ax        = ax,
    color     = 'lightblue',
    linewidth = 1,
    edgecolor = 'navy',
    label     = 'Original Boundaries'
)

plt.title('Distritos y Buffer')
plt.legend()
plt.axis('off')
(np.float64(-8272486.092246463),
 np.float64(-8197855.319316324),
 np.float64(4931921.070848358),
 np.float64(5006269.877978747))

Crear convex hulls:

gdf_3857['chull'] = gdf_3857.convex_hull

Visualización

Vamos a especificar parámetros para mejorar las visualizaciones:

plt.rcParams['figure.dpi'] = 150

Mapa temático:

## crear lienzo
fig, ax = plt.subplots(figsize = (10, 6))

## añadir mapa
gdf_3857.plot(
    column    = 'area',
    ax        = ax,
    legend    = True,
    cmap      = 'YlOrRd',
    edgecolor = 'black',
    linewidth = .5
)

## añadir elementos al mapa
plt.title('Distritos de NYC por área', fontsize = 16, fontweight = 'bold')
plt.axis('off')
plt.tight_layout()
plt.show()

Mapa interactivo con Folium:

gdf_3857.explore(
    column  = 'area',
    cmap    = 'YlOrRd',
    tooltip = ['area', 'dist_to_manhattan'],
    popup   = True,
    legend  = True
)
Make this Notebook Trusted to load map: File -> Trust Notebook

Relaciones espaciales

Las relaciones espaciales en Python se estrucutran: objeto1.relación(objeto2). Los métodos intersects, touches, crosses, contains, contains_properly… devuelven un Pandas Series de True/False.

Intersección con Manhattan:

## extraer Manhattan
manhattan_geom = gdf_3857.loc['Manhattan', 'geometry']

## intersección
intersects_manhattan = gdf_3857.intersects(manhattan_geom)

## filtrar los que intersecan
gdf_3857.loc[intersects_manhattan, :]
BoroCode Shape_Leng Shape_Area geometry area perim centroid boundary dist_to_manhattan buffer chull
BoroName
Queens 4 896344.047763 3.045213e+09 MULTIPOLYGON (((-8219461.955 4952778.645, -821... 4.928308e+08 360.377573 POINT (-8217436.761 4969318.623) MULTILINESTRING ((-8219461.955 4952778.645, -8... 19456.427458 POLYGON ((-8229469.036 4965425.107, -8229556.9... POLYGON ((-8231045.21 4944993.298, -8233479.86...
Brooklyn 3 741080.523166 1.937479e+09 MULTIPOLYGON (((-8222843.701 4950893.717, -822... 3.129684e+08 297.784594 POINT (-8231817.476 4960085.209) MULTILINESTRING ((-8222843.701 4950893.717, -8... 19586.763567 POLYGON ((-8245284.602 4956687.136, -8245296.1... POLYGON ((-8235790.052 4949053.259, -8237874.4...
Manhattan 1 359299.096471 6.364715e+08 MULTIPOLYGON (((-8238858.852 4965914.967, -823... 1.032201e+08 144.706085 POINT (-8233984.776 4979551.696) MULTILINESTRING ((-8238858.852 4965914.967, -8... 0.000000 POLYGON ((-8237068.209 4963506.169, -8237106.8... POLYGON ((-8240208.861 4965683.834, -8242627.6...
Bronx 2 464392.991824 1.186925e+09 MULTIPOLYGON (((-8226155.114 4982269.852, -822... 1.929250e+08 187.221157 POINT (-8222783.625 4990631.161) MULTILINESTRING ((-8226155.114 4982269.852, -8... 15755.010025 POLYGON ((-8233207.11 4988220.645, -8233189.61... POLYGON ((-8224095.48 4980733.074, -8225054.51...

Tocan a Manhattan:

gdf_3857.loc[gdf_3857.touches(manhattan_geom), :]
BoroCode Shape_Leng Shape_Area geometry area perim centroid boundary dist_to_manhattan buffer chull
BoroName
Queens 4 896344.047763 3.045213e+09 MULTIPOLYGON (((-8219461.955 4952778.645, -821... 4.928308e+08 360.377573 POINT (-8217436.761 4969318.623) MULTILINESTRING ((-8219461.955 4952778.645, -8... 19456.427458 POLYGON ((-8229469.036 4965425.107, -8229556.9... POLYGON ((-8231045.21 4944993.298, -8233479.86...
Brooklyn 3 741080.523166 1.937479e+09 MULTIPOLYGON (((-8222843.701 4950893.717, -822... 3.129684e+08 297.784594 POINT (-8231817.476 4960085.209) MULTILINESTRING ((-8222843.701 4950893.717, -8... 19586.763567 POLYGON ((-8245284.602 4956687.136, -8245296.1... POLYGON ((-8235790.052 4949053.259, -8237874.4...
Bronx 2 464392.991824 1.186925e+09 MULTIPOLYGON (((-8226155.114 4982269.852, -822... 1.929250e+08 187.221157 POINT (-8222783.625 4990631.161) MULTILINESTRING ((-8226155.114 4982269.852, -8... 15755.010025 POLYGON ((-8233207.11 4988220.645, -8233189.61... POLYGON ((-8224095.48 4980733.074, -8225054.51...

Comprobar que los centroides están dentro de su polígono:

gdf_3857['centroid'].within(gdf_3857.geometry)
BoroName
Staten Island    True
Queens           True
Brooklyn         True
Manhattan        True
Bronx            True
dtype: bool

Relaciones binarias (dos geometrías)

gdf_3857['intersection'] = gdf_3857.intersection(manhattan_geom)
gdf_3857['difference'] = gdf_3857.difference(manhattan_geom)

gdf_3857['intersection'].explore()
Make this Notebook Trusted to load map: File -> Trust Notebook

Verificaciones

gdf.is_empty
gdf.is_valid
gdf.is_simple
0    True
1    True
2    True
3    True
dtype: bool